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Hot Jupiters are among the best-studied exoplanets, but it is still poorly 
understood how their chemical composition and cloud properties vary 
with longitude. Theoretical models predict that clouds may condense 
onthe nightside and that molecular abundances can be driven out of 


equilibrium by zonal winds. Here we report a phase-resolved emission 
spectrum of the hot Jupiter WASP-43b measured from 5 um to 12 pm with 
theJWST's Mid-Infrared Instrument. The spectra reveal a large day-night 
temperature contrast (with average brightness temperatures of 1,524 + 35 K 
and 863 + 23 K, respectively) and evidence for water absorption at all orbital 
phases. Comparisons with three-dimensional atmospheric models show 
that both the phase-curve shape and emission spectra strongly suggest 

the presence of nightside clouds that become optically thick to thermal 
emission at pressures greater than -100 mbar. The dayside is consistent 
with a cloudless atmosphere above the mid-infrared photosphere. 
Contrary to expectations from equilibrium chemistry but consistent with 
disequilibrium kinetics models, methane is not detected on the nightside 
(2c upper limit of 1-6 ppm, depending on model assumptions). Our results 
provide strong evidence that the atmosphere of WASP-43b is shaped by 
disequilibrium processes and provide new insights into the properties of the 
planet's nightside clouds. However, the remaining discrepancies between 
our observations and our predictive atmospheric models emphasize the 
importance of further exploring the effects of clouds and disequilibrium 
chemistry in numerical models. 


Hot Jupiters are tidally synchronized to their host stars, with vast 
differences in irradiation between the dayside and nightside. Previ- 
ous observations with the Hubble Space Telescope (HST) and the 
Spitzer Space Telescope show that these planets have cooler night- 
sides and weaker hotspot offsets than expected from cloud-free 
three-dimensional models'^. The main mechanism believed to be 
responsible for this behaviour is the presence of nightside clouds, 
which would hide the thermal flux of the planet and lead to a sharp lon- 
gitudinal gradient in brightness temperature^^*?, Other mechanisms 
have been proposed, such as the presence of atmospheric drag due 
to hydrodynamic instabilities or magnetic coupling!" ©, super-stellar 
atmospheric metallicity”, or interaction between the deep winds 


and the photosphere!5, but these mechanisms are less universal than 
the cloud hypothesis", 

WASP-43b, a hot Jupiter with an orbital period ofjust 19.5 h (ref. 19), 
isanidealtarget for thermal phase-curve observations. Its host star is 
aK7 main-sequence star 87 pc away with metallicity closeto solar and 
weak variability”°. Previous measurements ofthe planet's orbital phase 
curve in the near-infrared have revealed a large temperature contrast 
between the dayside and nightside hemispheres, broadly consistent 
with the presence of nightside clouds*”’”, which could be composed of 
magnesium silicates (Mg,SiO,/MgSiO,) and other minerals (for exam- 
ple, MnS, Na,S, metal oxides)”**. Owing to the low nightside flux, the 
exact temperature and cloud properties were challenging to determine 
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Fig. 1| A visualization of the observed light curves and the resulting emission 
spectra. a, The observed spectroscopic light curves binned to a 0.5 um 
wavelength resolution and after systematic noise removal, following the Eureka! 
vi methods. The first 779 integrations have been removed from this figure and 
our fits as they were impacted by strongly decreasing flux. Wavelengths longer 
than 10.5 pm marked with a hatched region were affected by the ‘shadowed 
region effect’ (Methods) and could not be reliably reduced. b, The observed 
band-integrated light curve after systematic noise removal (grey points) and 
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binned data with a cadence of 15 min (black points, with error bars smaller than 
the point sizes), compared with the best-fitting astrophysical model (red line). 
c,d, The measured dayside (c) and nightside (d) emission spectra are shown with 
black points and 1c error bars, and black-body curves (dotted line denoted as ‘BB’, 
assuming a PHOENIX” 75 model for the star) are shown to emphasize planetary 
spectral features with black-body temperatures estimated by eye to match the 
continuum flux levels. Wavelengths longer than 10.5 jum were affected by the 
shadowed region effect and are unreliable. 


from previous observations*?^?*, With the mid-infrared capabilities of 
theJWST, we have the opportunity to measure the phase-resolved ther- 
mal spectrum with unprecedented sensitivity, particularly onthe cold 
nightside. We observed a full orbit of WASP-43b in the 5-12 um range 
with the JWST’s Mid-Infrared Instrument (MIRI)” in low-resolution 
spectroscopy (LRS)? slitless mode on 1 and 2 December 2022, as part 
of the Transiting Exoplanet Community Early Release Science Pro- 
gram (JWST-ERS-1366). This continuous observation lasted 26.5 hata 
cadence of 10.34 s (9,216 integrations) and included a full phase curve 
with one transit and two eclipses. 


Results 

We performed multiple independent reductions and fits to these 
observations (see ‘Data reduction pipelines’ and ‘Light-curve fitting’ 
in Methods) to ensure robust conclusions. Our analyses all identified 
a strong systematic noise feature from 10.6 um to 11.8 um, the source 
of which is still unclear, and we were unable to adequately detrend 
these 10.6-11.8 jum data (see 'Shadowed region effect' in Methods). 
As shown in Extended Data Fig. 1, we also found that larger wavelength 
bins were required to accurately estimate our final spectral uncertain- 
ties (see ‘Spectral binning’ in Methods). As a result, our final analyses 
consider only the 5-10.5 jum data, which we split into 11 channels with 
a constant 0.5 uim wavelength spacing. Similar to the MIRI commis- 
sioning time-series observations, our data show a strong downwards 
exponential ramp in the first -60 min and a weaker ramp throughout 
theobservation? (Extended Data Fig. 2). To minimize correlations with 
the phase variations, we removed the initial strong ramp by excluding 
thefirst 779 integrations (134.2 min) and then fitted a single exponential 
ramp modelto the remaining data. A single ramp effectively removed 


the systematic noise, with the broadband light curve showing scat- 
ter -1.25x the expected photon noise, while the spectroscopic light 
curves reach aslow as -1.1x the photon limit, probably due to improved 
decorrelation of wavelength-dependent systematics. Figure 1 shows 
the spectral light curves, broadband light curve, dayside spectra and 
nightside spectra from our fiducial reduction and fit. 

From our Eureka! v1 analysis (Methods), we measure a broad- 
band (5-10.5 um) peak-to-trough phase variation of 4,180 + 33 ppm 
with an eclipse depth of 5,752 +19 ppm and a nightside flux of 
1,636 + 37 ppm. Assuming a PHOENIX stellar model and marginalizing 
over the published stellar and system parameters”, the broadband 
dayside brightness temperature is 1,524 + 35 K while the nightside is 
863 + 23 K. This corresponds to a day-night brightness temperature 
contrast of 659 +19 K, in agreement with the large contrasts previ- 
ously observed*?'???, The phase variations are well fitted by a sum of 
two sinusoids (the first and second harmonics), with two sinusoids 
preferred over a single sinusoid at 160 (see ‘Determining the number 
of sinusoid harmonics’ in Methods) for the broadband light curve. The 
peak brightness of the broadband phase curve occurs at 7.34 + 0.38? E 
from the substellar point (although individual reductions find offsets 
ranging from 7.34° E to 9.60° E), while previous studies have found 
offsets of 12.3 + 1.0? E for HST Wide Field Camera 3’s (WFC3) 1.1-1.7 um 
bandpass”, offsets ranging from 4.4? E to 12.2? E for Spitzer Infra- 
Red Array Camera’s (IRAC) 3.6 pm filter??? and offsets ranging from 
10.4? E to 21.1? E for Spitzer/IRAC's 4.5 jum filter*””>°*"*, Overall, 
these broadband data represent roughly an order of magnitude in 
improved precision onthe eclipse depth (6x), phase-curve amplitude 
(6x) and phase-curve offset (10x) over individual Spitzer/IRAC 4.5 ym 
observations of the system??*?" thisimprovementis largely driven by 
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Fig. 2 | A comparison of the observed 5-10.5 pm light curve with GCM 
simulations. The black points show the temporally binned broadband light 
curve. The solid lines represent modelled phase curves derived from the 

31GCM simulations, integrated over the same wavelength range as the data, and 
separated into two groups based on the inclusion of clouds. The cloudless GCMs 
(red lines) simulated completely cloud-free skies, whereas the cloudy GCMs 
(blue lines) included at least some clouds on the nightside of the planet. The red 
and blue shaded areas span the range ofall the cloudless and cloudy simulations, 
respectively, with the spread of values owing to differences in the various model 
assumptions and parameterizations. On average, the cloudless GCM phase 
curves have a maximum planet-to-star flux ratio of 5,703 ppm and a minimum 

of 2,681 ppm. This matches the observed maximum of the phase curve well but 
does not match its observed minimum at 1,636 * 37 ppm. On average, the cloudy 
GCM phase curves have a maximum of 5,866 ppm and a minimum of 1,201 ppm, 
in better agreement with the observed nightside emission, but their spread 

of maximum values is much larger than the cloudless simulations. The cloudy 
models are able to suppress the nightside emission and better match the data; 
however, not all cloud models fit equally well and those with the optically thickest 
nightside clouds suppress too much emission. The models do not include the 
eclipse signals (phases -0.5 and 0.5) or transit signal (phase 0.0). 


theJWST'slarger mirror (45x), about 12x less pointing jitter (per axis), 
about 4x improved stability in the width ofthe point spread function 
(PSF) along each axis and MIRI's much broader bandpass. 


Model interpretation 

To interpret the measurements, we compared the observations with 
synthetic phase curves and emission spectra derived from general 
circulation models (GCMs). Simulations were gathered from five dif- 
ferent modelling groups, amounting to 31 separate GCM realizations 
exploring a range of approaches and assumptions. Notably, in addi- 
tion to cloud-free simulations, the majority of the GCMs modelled 
clouds with spatial distributions that were either fully predicted^?9? 
or simply limited to the planet's nightside*. For the predictive cloud 
models, simulations favoured warmer, clearer daysides with cooler, 
cloudier nightsides, butthe precise distributions varied with assump- 
tions regarding cloud physics and compositions. In general, models 
with smaller cloud particles or extended vertical distributions tended 
to producethicker clouds at the pressures sensed by the observations. 
Details of the different models are provided in Methods. 

Despite fundamental differences in the models and the param- 
eterizations they employ, simulated phase curves derived from models 
that include cloud opacity on the planet's nightside provide a better 
match to the observed nightside flux compared with the clear simula- 
tions (Fig. 2). In contrast, the observed dayside fluxes (180? orbital 
phase) were matched similarly well by models with and without clouds. 
This implies the presence of widespread clouds preferentially on the 
planet's nightside with cloud optical thicknesses sufficientto suppress 
thermal emission and cool the thermal photosphere. Specifically, 
models with integrated mid-infrared cloud opacities of roughly 2-4 
abovethe300 mbar level (that is, blocking -87-98% of the underlying 
emission), best match the observed nightside flux. 


Including nightside clouds also improved the agreement with 
the measured hotspot offset (7.34+ 0.38° E). While cloudless models 
all produced eastward offsets greater than 16.6? (25.5? on average), 
simulations with clouds had offsets as low as 7.6° (with a mean of 16.4°). 
These reduced offsets were associated with decreases in the eastwards 
jet speeds of up to several kilometres per second, with maximum winds 
of roughly 2.0-2.5 km s ! providing the best match (see Extended Data 
Table 1 for further details). This modelled jet-speed reduction is prob- 
ably due toa disruption in the equatorwards momentum transport** 
brought about by nightside clouds^?^?*, However, the resulting range 
of offsets seen in the suite of models suggests that this mechanism is 
quite sensitive to the details of cloud models, and other modelling 
factors (for example, atmospheric drag!""?!5 radiative timescales?) 
probably still play an important role. 

A comparison of the observed and modelled emission spectra 
further suggests that the majority ofthe cloud thermal opacity must 
be confined to pressures greater than -10-100 mbar, because the 
presence of substantial cloud opacity at lower pressures dampens 
the modelled spectral signature amplitude below what is observed 
(Fig. 3). No distinct spectral signatures indicative of the cloud com- 
position were evident in the observations. While no single GCM can 
match the emission spectra at all phases, spectra corresponding to 
nightside, morning and evening terminators appear qualitatively 
similar to GCM results that are intermediate between clear and 
cloudy simulations. In contrast, the absorption features indicative 
of water vapour (between -5 um and 8.5 jum) seen in the dayside 
emission spectrum are more consistent with an absence of cloud 
opacity atthese mid-infrared wavelengths. Altogether, these findings 
represent new constraints on the spatial distribution and opacity of 
WASP-43b's clouds. 

Wefurther characterized the chemical composition of WASP-43b's 
atmosphere by applying a suite of atmospheric retrieval frameworks to 
the phase-resolved emission spectra. The retrievals spanned a broad 
range of model assumptions, including free chemical abundances 
versus equilibrium chemistry, different temperature profile param- 
eterizations and different cloud models (see 'Atmospheric retrieval 
models' in Methods). Despitethese differences, the retrievals yielded 
consistent results for both the chemical and thermal constraints. We 
detected water vapour across the dayside, nightside, morning and 
evening hemispheres, with detection significances of up to ~3-40 
(Extended Data Fig. 3 and Extended Data Tables 2 and 3). The retrieved 
abundances of HO largely lie in the 10-10? ppm range for all four phases 
and for all the retrieval frameworks (Fig. 4 and Extended Data Fig. 4), 
broadly consistent with the value expected for a solar composition 
(500 ppm) as wellas previous observations”. 

We also searched for signatures of disequilibrium chemistry in 
the atmosphere of WASP-43b. While CH, is expected to be present on 
the nightside under thermochemical equilibrium conditions, we did 
not detect CH, at any phase (Fig. 4). In the pressure range probed by 
the nightside spectrum (1-10? bar; Extended Data Fig. 5), the equi- 
librium abundance of CH, is expected to vary between -1 ppm and 
100 ppm for a solar C/O ratio?, compared with our 9596 upper limits 
of 1-6 ppm (Extended Data Table 2). The upper limits we place on the 
nightside CH, abundance are more consistent with disequilibrium 
models that account for vertical and horizontal transport^?^*^5, In 
particular, two-dimensional photochemical models and GCMs pre- 
dict the strongest depletion of CH, on the nightside due to strong 
zonal winds (>1 km s?) transporting gas-phase constituents around the 
planet faster than the chemical reactions can maintain thermochemical 
equilibrium, thus ‘quenching’ and homogenizing the global compo- 
sition at values more representative of dayside conditions (see also 
refs. 39-42). We note, however, that alow atmospheric C/O ratio and/ 
or clouds at photospheric pressures could also lead to a non-detection 
of CH,. We also searched for signatures of NH}, which is predicted to 
havea volume mixing ratio less than 0.1-1 ppm in both equilibrium and 
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Fig. 3 | A comparison of the observed and modelled spectra at different 
phases. a-c, The observed emission spectrum with lo error bars at phases 0.0 
(a), 0.25 (b), 0.5 (c) and 0.75 (d), along with select modelled spectra derived 
from different cloudy and cloudless GCMs (described in Methods and listed 
in Extended Data Table 1). Although absolute brightness temperatures differ 
appreciably between models owing to various GCM assumptions, differences 
inthe relative shape ofthe spectra are strongly dependent on the cloud and 
temperature structure found in the GCMs (Extended DataFig. 7). Models 

with more isothermal profiles (like RM-GCM) or thick clouds at pressures of 
<S10-100 mbar (like THOR cloudy, Generic PCM with 0.1 um cloud particles) 
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produce flatter spectra, while clearer skies yield stronger absorption features. 
The observed spectra from the nightside and terminators appear muted 
compared with the clear-model spectra, suggesting the presence of at least some 
clouds or weak vertical temperature gradients at pressures of 10-100 mbar. In 
contrast, the spectral structure produced by water vapour opacity (indicated by 
the purple shading) appears more consistent with models lacking clouds at these 
low pressures on the dayside. Under equilibrium chemistry, methane would also 
show an absorption feature at -7.5-8.5 um (shaded pink) for the colder models at 
phases 0.0 and 0.75. Finally, the median retrieved spectrum and lo contours from 
the HyDRA retrieval are shown in grey. 


disequilibrium chemistry models, and find that the results are incon- 
clusive and model-dependent with the current retrieval frameworks. 

Given the strong evidence for clouds from comparison with GCMs, 
we also searched for signatures of clouds inthe atmospheric retrieval. 
Formally, the retrievals do not detect clouds with statistical signifi- 
cance, indicating that strong spectral features uniquely attributable 
to condensates are not visible in the data (see Atmospheric retrieval 
models' in Methods and Extended DataFig. 6). However, the retrievals 
may mimic the effects of cloud opacity with a more isothermal tem- 
perature profile, as both tend to decrease the amplitude of spectral 
features, but the cloud-free, more isothermal temperature profile 
requires fewer free parameters and is therefore statistically favoured. 
Indeed, while the retrieved temperature profiles on the dayside and 
evening hemispheres agree well with the hemispherically averaged 
temperature profiles from the GCMs, they are more isothermal than 
the GCM predictions for the nightside and morning hemispheres 


(Extended Data Fig. 7). This discrepancy may hint at the presence of 
clouds on the nightside and morning hemispheres, consistent with 
thelocations of clouds found in the GCMs. 

Taken together, our results highlight the unique capabilities of 
JWST/MIRI for exoplanet atmosphere characterization. Combined with 
arange of atmospheric models, the observed phase curve and emission 
spectra provide strong evidence that the atmospheric chemistry of 
WASP-43b is shaped by complex disequilibrium processes and provide 
new constraints on the optical thickness and pressure of nightside 
clouds. However, while cloudy GCM predictions match the data better 
than cloud-free models, none of the simulations simultaneously repro- 
duced the observed phase curve and spectra within measured uncer- 
tainties. These remaining discrepancies underscore the importance of 
further exploring the effects of clouds and disequilibrium chemistry in 
numerical models, as JWST continues to place unprecedented obser- 
vational constraints on smaller and cooler planets. 
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Fig. 4 | Atmospheric spectral retrievals for frameworks with free chemistry. 
a, Temperature profile contours (68% confidence) constrained by the retrievals 
at each orbital phase (see legends). All frameworks produced consistent non- 
inverted thermal profiles that are consistent with two-dimensional radiative- 
convective equilibrium and photochemical models along the equator” (black 
curves) over the range of pressures probed by the observations (black bars). 

b, H;O abundance posterior distributions (volume mixing ratios). The shaded 
areas denote the span of the 68% confidence intervals. The green and blue bars on 


3,000 
logXc, 


each panel denote the abundances predicted by equilibrium and disequilibrium 
chemistry solar-abundance models”, respectively, at the pressures probed 

by the observations (1-10? bar, approximately). c, The same as in b but for 

CH,. Theretrieved water abundances are consistent with either equilibrium or 
disequilibrium chemistry estimations for solar composition (500 ppm), whereas 
the retrieved upper limits to the CH, abundance are more consistent with 
disequilibrium chemistry predictions. 


Methods 

Observations and quality of the data 

We observed a full orbit of WASP-43b with the JWST MIRI LRS slitless 
mode as a part of JWST-ERS-1366. We performed target acquisition 
with the F1500W filter and used the SLITLESSPRISM subarray for the 
science observation. The science observation was taken between 1 
December 2022 at 00:54:30 UT and 2 December 2022 at 03:23:36 UT, 
for a total of 26.5 h. We acquired 9,216 integrations, which were split 
into3 exposures and 10 segments per exposure. Each integration lasts 
10.34 s and is composed of 64 groups, with 1 frame per group. The 
LRS slitless mode reads an array of 416 x 72 pixels on the detector (the 
SLITLESSPRISM subarray) and uses the FASTRI readout mode, which 
introduces an additional reset between integrations. 

Owing to the long duration of the observation, two high-gain 
antenna moves occurred 8.828 h and 17.661 h after the start of the 
science observation. They affect only a couple of integrations that we 
removed from the light curves. A cross-shaped artefact is present on 
thetwo-dimensional images at the short-wavelength end dueto light 
scattered by detector pixels”. It is stable over the duration ofthe obser- 
vationbutit contaminates the background and the spectraltrace up to 
-6 um. This ‘cruciform’ artefact is observed in all MIRI LRS observations; 
a dedicated analysis is underway to estimate and mitigate its impact. 

In the broadband light curve, the flux decays by ~0.1% during the 
first 60 min and continues to decay throughout the observation. This 
ramp is well modelled with 1 or 2 exponential functions after trimming 
the initial -780 integrations. Without trimming any data, at least two 
ramps are needed. In addition, a downwards linear trend in flux is 
observed overthe whole observation with a slope of -39 ppm per hour. 
Thesetwo types of drift also appear in the spectroscopic light curves. 
Theexponential ramp amplitude in thefirst 60 min changes with wave- 
length from -0.67% in the 5-5.5 um bin (downwards ramp) to +0.26% 
inthe10-10.5 um bin (upwards ramp). The ramp becomes upwards at 
wavelengths longer than 7.5 jum and its timescale increases to more 
than 1h at wavelengths longer than 10.5 um. The slopes as a function 
of wavelength vary from -16 ppm to -52 ppm, all downwards. Such 
drifts (initial ramp and linear or polynomial trend) are also observed in 
other MIRILRS time-series observations? butthe strength ofthe trends 
differ for each observation. In these WASP-43b observations, we note 
that their characteristic parameters vary smoothly with wavelength, 
which may help identify their cause and build correction functions. 


Over the course of the observation, the position of the spectral 
trace onthe detector varies by 0.0036 pixels RMS (0.027 pixels peak to 
peak) in the spatial direction, and the Gaussian standard deviation of 
the spatial PSF varies by 0.00069 pixels RMS (root mean square; 0.0084 
pixels peak to peak) following a sharp increase by 0.022 pixels during 
the first 600 integrations. Depending on the wavelength bin, that spatial 
drift causes noise at the level of 7-156 ppm, while variations in the PSF 
width cause noise at the level of 4-54 ppm (these numbers are obtained 
froma linear decorrelation). Overall, the MIRI instrument used in LRS 
slitless moderemains remarkably stable over this 26.5-h-long continu- 
ous observation and the data are of exquisite quality. 

The noise in the light curve increases sharply at wavelengths 
beyond 10.5 um and the transit depths obtained at these long wave- 
lengths by different reduction pipelines are discrepant. These wave- 
lengths were not used in the retrieval analyses and the final broadband 
light curve. The causeisunknownbutit might be related to the fact that 
this region of the detector receives different illumination before the 
observation“ (see ‘Shadowed region effect’ below for more details). 


Data reduction pipelines 

Eureka! v1 reduction. The Eureka! v1 reduction made use of version 
0.9 of the Eureka! pipeline“, CRDS version 11.16.16 and context 1018, 
and jwst package version 1.8.3*°. The gain value of 5.5 electrons per 
data number obtained from these CRDS reference files is known to be 
incorrect, and the actual gain is estimated to be -3.1 electrons per data 
number although the gain may be wavelength dependent (S. Kendrew, 
private communication). A new reference file reflecting the updated 
gain is under development at STScl, which will improve the accuracy 
of photon-noise calculations. For the rest of this analysis, we assume 
aconstant gain of 3.1 electrons per data number. The Eureka! control 
files and Eureka! parameter files files used in these analyses are avail- 
able for download (https://doi.org/10.5281/zenodo.10525170) and are 
summarized below. 

Eureka! makes use of the jwst pipeline for stages 1 and 2, and 
both stages were run with their default settings, with the exception 
of increasing the stage 1 jump step’s rejection threshold to 8.0 and 
skipping the photom step in stage 2 because it is not necessary and 
can introduce additional noise for relative time-series observations. 
In stage 3 of Eureka!, we then rotated the MIRI/LRS slitless spectra 90° 
anticlockwise so that wavelength increases from left to right like the 
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other JWST instruments to allow for easier reuse of Eureka! functions. 
We then extracted pixels 11-61 in the new y direction (the spatial direc- 
tion) and 140-393 in the new x direction (spectral direction); pixels 
outside of these ranges primarily contain noise that is not useful for our 
reduction. Pixels marked as DO NOT USE' in the DQ array were then 
masked as were any other unflagged NaN or inf pixels. A centroid was 
then fit to each integration by summing along the spectral direction 
andfittingthe resulting one-dimensional profile with a Gaussian func- 
tion; the centroid from the first integration was used for determining 
aperture locations, while the centroids and PSF widths from all integra- 
tions were saved to be used as covariates when fitting the observations. 

Our background subtraction method is tailored to mitigate sev- 
eral systematic effects unique to the MIRI instrument. First, MIRI/ 
LRS observations exhibit a ‘cruciform artefact” at short wavelengths 
caused by scattered light within the optics; this causes bright rays of 
scattered light which must be sigma-clipped to avoid over-subtracting 
the background. In addition, MIRI/LRS observations show periodic 
noise in the background flux, which drifts with time” as well as 1/f 
noise*, which leads to correlated noise in the cross-dispersion direc- 
tion; as a result, background subtraction must be performed indepen- 
dently for each integration and column (row in MIRI's rotated reference 
frame). Furthermore, in both these observations and the dedicated 
background calibration observations from JWST-COM/MIRI-1053, we 
found that there was a linear trend in the background flux, with the 
background flux increasing with increasing row index (column index 
in MIRI's rotated reference frame). To robustly removethis feature, we 
found that it was important to either (1) use the mean from an equal 
number of pixels on either side of the spectral trace for each column 
andintegration, or (2) usealinear background modelfor each column 
and integration; we adopted the former as it resulted in less noisy light 
curves. To summarize, for each column in each integration we sub- 
tracted the mean of the pixels separated by 211 pixels from the centre 
ofthe spectral trace after first masking So outliers in that column. 

To compute the spatial profile for the optimal extraction of the 
source flux, we calculated a median frame, sigma-clipping 5o outliers 
along the time axis and smoothing along the spectral direction using a 
7-pixel-wide boxcar filter. During optimal extraction, we only used the 
pixels within 5 pixels ofthe fitted centroid and masked pixels that were 
100 discrepant with the spatial profile. Background exclusion regions 
ranging from 9 to 13 pixels and source aperture regions ranging from 
4 to 6 pixels were considered, but our values of 11 and 5 pixels were 
selected as they produced the lowest median absolute deviation light 
curves before fitting. 


Eureka! v2 reduction. The Eureka! v2 reduction followed the same pro- 
cedureasthe Eureka! v1 reduction except for the following differences. 
First, this reduction made use of version 1.8.1 of thejwst pipeline. For 
stage 1, we instead used a cosmic ray detection threshold of 5and used a 
uniform ramp fitting weighting. For stage 2, we performed background 
subtraction using columns away from the trace on the left and on the 
right and subtracted the background for each integration”. Stage 3 
was identical to Eureka! v1 reduction. 


TEATRO reduction. We processed the data using the Transiting Exo- 
planet Atmosphere Tool for Reduction of Observations (TEATRO) that 
runs thejwst package, extracts and cleans the stellar spectra and light 
curves, and runs light-curve fits. We used thejwst package version 1.8.4, 
CRDS version 11.16.14 and context 1019. We started from the 'uncal' files 
and ran stages 1 and 2 of the pipeline. For stage 1, we set a jump rejec- 
tionthreshold of 6, turned offthe jump.flag 4 neighbors' parameter 
and used the default values for all other parameters. For stage 2, we 
ran only the ‘AssignWcsStep, ‘FlatFieldStep’ and ‘SourceTypeStep’; no 
photometric calibration was applied. The next steps were made using 
our ownroutines. We computed the background using two rectangular 
regions, one on each side of the spectral trace, between pixels 13 and 


27 and between pixels 53 and 72 in the spatial direction, respectively. 
We computed the background value for each row (rows are along the 
spatial direction) in each region using a biweight location, averaged the 
two values and subtracted it from the full row. This background sub- 
traction was done for each integration. Then, we extracted the stellar 
flux using aperture photometry by averaging pixels between 33 and 42 
ineachrowto obtain the stellar spectrum at each integration. We also 
averaged pixels between 33 and 42 in the spatial direction and between 
5 um and 10.5 pm in the spectral direction to obtain the broadband 
flux. We averaged the spectra in 11 0.5-um-wide wavelength channels. 
For each channel and for the broadband light curve, we normalized 
the light curve using the second eclipse as a reference flux, computed 
a running median filter using a 100-point window size, and rejected 
points that were more than 3g away from that median using a 5-iteration 
sigma-clipping. To limit the impact of the initial ramp on the fitting, 
we trim the first 779 integrations from the broadband light curve and 
a similar number of integrations for each channel (the exact number 
depends onthe channel). Finally, we subtracted 1 from the normalized 
light curves to have the secondary eclipse flux centred on O. These 
cleaned light curves were used for phase curve, eclipse and transit fits. 


SPARTA reduction. We reduced the data with the open-source Simple 
Planetary Atmosphere Reduction Tool for Anyone (SPARTA), first 
introduced in ref. 48 to analyse the MIRI LRS phase curve of GJ 1214b. 
We started from the uncalibrated data and proceeded all the way to 
the final results without using any code from the jwst or Eureka! pipe- 
lines. In stage 1, we started by discarding the first five groups as well 
as the last group, because these groups show anomalies due to the 
reset switch charge decay and the last-frame effect. We fitted a slope 
to the up-the-ramp reads in every pixel of every integration in every 
exposure. We calculated the residuals of these linear fits, and for every 
pixel, we computed a median residual for every group across all inte- 
grations. This ‘median residual’ array has dimensions Na, X Nrows X Neots: 
This array was subtracted from the original uncalibrated data and the 
up-the-ramp fit was redone, this time without discarding any groups 
except those that were more than 5e away from the best-fit line. Such 
outliers, which may be due to cosmic rays, were discarded and the fit 
recomputed until convergence. This procedure straightens out any 
nonlinearity inthe up-the-ramp reads that is consistent across integra- 
tions, such as the reset switch charge decay, the last-frame effect or 
inaccuracies in the nonlinearity coefficients. After up-the-ramp fitting, 
weremoved the background by removing the mean of columns 10-24 
and 47-61 (inclusive, zero-indexed) for every row of every integration. 
As these two regions are of equal size and equally distant from the 
trace, any linear spatial trend in the background is naturally removed. 

In the next step, we computed pixel-wise median image over all 
integrations. This median image was used asa template to determine 
the position of the trace in each integration, by shifting and scaling 
the template until it matched the integration (and minimizes the x°). 
It was also used as the point spread profile for optimal extraction, 
after shifting in the spatial direction by the amount calculated in the 
previous step. Outliers more than 5o discrepant from the model image 
(which may be cosmic rays) were masked, and the optimal extraction 
was repeated until convergence. The z-scores image (image minus 
model image all divided by expected error, including photon noise 
and read noise) have a typical standard deviation of 0.88, compared 
with a theoretical minimum value of 1, indicating that the errors are 
being overestimated. 

After optimal extraction, we gathered all the spectra and posi- 
tions into one file. To reject outliers, we created a broadband light 
curve, detrended it by subtracting a median filter with a width 100 
times less than the total data length and rejected integrations greater 
than 4g away from 0 (which may be cosmic rays). Sometimes only cer- 
tain wavelengths of an integration are bad, not the entire integration. 
We repaired these by detrending the light curve at each wavelength, 
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identifying 40 outliers and replacing them with the average of their 
twoimmediate temporal neighbours. 


Spectral binning 

To investigate the effects of spectral binning, we utilized the MIRI 
time-series commissioning observations of the transit of L168-9b 
(JWST-COM/MIRI-1033; ref. 29). L168-9b was chosen to have a clear 
transit signal while also having no detectable atmospheric signatures 
expected in its mid-infrared transmission spectrum; as a result, the 
observed scatter in the transmission spectrum can be used as an inde- 
pendent measurement of the uncertainties in the transit depths. The 
same procedure cannot be done on our WASP-43b science observations 
as there may be detectable atmospheric signatures. 

Following the Eureka! reduction methods described by ref. 29, 
we tried binning the L168-9b spectroscopic light curves at different 
resolutions and compared the observed standard deviation of the 
transmission spectrum with the median of the transit depth uncer- 
tainties estimated from fitting the spectral light curves. As shown in 
Extended Data Fig. 1, the uncertainties in the pixel-level light curves 
underestimate the scatter in the transmission spectrum by a factor of 
about two. Because pairs of rows (in MIRI's rotated reference frame) 
are reset together, it is reasonable to assume that there could be 
odd-even effects that would average out if combining pairs of rows; 
indeed, there do appear to be differences in the amplitude of the ini- 
tial exponential ramp feature between odd and even rows. However, 
combining pairs of rows still leads to appreciable underestimation of 
thescatter inthetransmission spectrum. Interestingly, the underesti- 
mation ofthe uncertainties appears to decrease with decreasing wave- 
length resolution. This is likely explained by wavelength-correlated 
noise that gets averaged out with coarse binning. A likely culprit for 
this wavelength-correlated noise may be the 390 Hz periodic noise 
observed in several MIRI subarrays, which causes clearly structured 
noise with a period of -9 rows? (M. Ressler, private communication); 
this noise source is believed to be caused by MIRI's electronics and pos- 
sible mitigation strategies are still under investigation. Until the source 
of the excess wavelength-correlated noise is definitively determined 
andanoise mitigation methodis developed, we recommend that MIRI/ 
LRS observations should be binned to a fairly coarse spectral resolution 
as this gives better estimates of the uncertainties and also gives spectra 
that are closer to the photon-limited noise regime. However, we caution 
against quantitative extrapolations of the uncertainty underestimation 
to other datasets; because we do not know the source of the excess 
noise, we do not know howit might change with different parameters 
such as groups per integration or stellar magnitude. 

Ultimately, for each reduction method, we binned the spectra 
downto aconstant 0.50-pm-wavelength grid spanning 5-12 um, giving 
atotal of 14 spectral channels. However, as is described below, we only 
end up using the 11 spectral channels spanning 5-10.5 um for science. 
This 0.5-um-binning scheme combines 7 wavelengths for the short- 
est bin and 25 wavelengths for the longest bin, which has the added 
benefit of binning down the noise at longer wavelengths where there 
are fewer photons. However, even for this coarse of abinning scheme, 
we do expectthere to be some additional noise beyond our estimated 
uncertainties on the spectrum of WASP-43b (Extended Data Fig. 1). 
Asthe structure of this noise source is not well understood nor is the 
extent to which our error bars are underestimated, our best course of 
action wasto consider error inflation when performing spectroscopic 
inferences (described in more detail below). 


Light-curve fitting 

Detrending the initial exponential ramp. As with other MIRI/LRS 
observations”, our spectroscopic light curves showed a strong expo- 
nential ramp at the start of the observations. As shown in Extended 
Data Fig. 2, the strength and sign ofthe ramp varies with wavelength, 
changing froma strong downwards ramp at 5 um to a nearly flat trend 


around 8 pm, and then becoming an upwards ramp towards longer 
wavelengths. From 10.6 um to 11.8 um, the ramp timescale became 
muchlonger and the amplitude ofthe ramp became much stronger; this 
region of the detector was previously discussed" and is discussed in 
more detail below. In general, most ofthe ramp's strength had decayed 
within -30-60 min, butatthe precision of our data, the residual ramp 
signature still had an important impact on our nightside flux meas- 
urements due to the similarity of the ramp timescale with the orbital 
period. Unlike inthe case of the MIRI/LRS commissioning observations 
of L168-9b”’, we were not able to safely fit the entire dataset with a small 
number of exponential ramps. When fitting the entire dataset, we found 
that non-trivial choices about the priors for the ramp amplitudes and 
timescales resulted in significantly different spectra at phases 0.75 
(morning hemisphere) and especially 0.0 (nightside); because the 
dayside spectrum is measured again near the end of the observations, 
it was less affected by this systematic noise. 

Ultimately, we decided to conservatively discard the first 779 
integrations (134.2 min), leaving only one transit duration of baseline 
before the first eclipse ingress began. After removing the initial 779 
integrations, we found that a single exponential ramp model with 
broad priors that varied freely with wavelength was adequate to remove 
the signature. In particular, after removing the first 779 integrations 
we found that our dayside and nightside emission spectra were not 
significantly affected by (1) fitting two exponential ramps instead of 
one, (2) adjusting our priors onthe ramp timescale to exclude rapidly 
decaying ramps with timescales >15 d! instead of >100 d”, (3) putting 
auniform prior onthe inverse timescale instead of the timescale, or (4) 
altering the functional form of the ramp by fitting for an exponential 
to which the time was raised. After removing the first ~2 h, we also 
found that the ramp amplitude and timescale did not vary strongly 
with wavelength (excluding the ‘shadowed region’ described below), 
although fixing these parameters to those fitted to the broadband light 
curve affected several points in the nightside spectrum by more than 
1o; we ultimately decided to leave the timescale and amplitude to vary 
freely with wavelength as there is no a priori reason to assume that they 
should be equal. With careful crafting of priors, it appeared possible 
to get results similar to our final spectra while removing only the first 
few integrations, but trimming more integrations and only using a 
single exponential ramp model required fewer carefully tuned prior 
assumptions for which we have little physical motivation. 


Shadowed region effect. As was described in ref. 44, we also identi- 
fied a strong discontinuity in the spectroscopic light curves span- 
ning pixel rows 156-220 (10.6-11.8 um) in these observations. In this 
range, the temporal behaviour of the detector abruptly changes to a 
large-amplitude, long-timescale, upwards ramp that appears to slightly 
overshoot before decaying back down and approaching an equilibrium. 
These pixels coincide with a region of the detector between the Lyot 
coronagraph region and thefour-quadrant phase mask region, which 
is unilluminated except when the dispersive element is in the optical 
path; as a result, we have taken to calling this unusual behaviour as 
the ‘shadowed region effect’. Strangely, not all MIRI/LRS observations 
show this behaviour, with the MIRI/LRS commissioning time-series 
observations? and the GJ 1214b phase-curve observations ^? showing 
no such effect. In fact, we know of only two other observations that 
show similar behaviour: the observation of the transit of WASP-80b 
(JWST-GTO-1177; T. Bell, private communication) and the observa- 
tion of the phase curve of GJ 367b (JWST-GO-2508; M. Zhang, private 
communication). Informatively, the eclipse observation of WASP-80b 
taken 36 h after the WASP-80b transit using the same observing pro- 
cedure (JWST-GTO-1177; T. Bell, private communication) did not show 
the same shadowed region effect, indicating that the effect is unlikely 
to be caused by stray light from nearby stars or any other factors that 
stayed the same between those two observations. Our best guess at 
this point is that the effect is related to the illumination history of the 
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detector and the filter used by the previous MIRI observation (because 
the detector is illuminated at all times, even when it is not in use), but 
thisis still under investigation and at present there is no way of predict- 
ing whether or not an observation will be impacted by the shadowed 
region effect. It is important to note, however, that from our limited 
knowledge at present that the shadowed region effect appears to be 
either present or not, with observations either strongly affected or 
seemingly completely unaffected. 

Using the general methods described in the Eureka! v1 fit, we 
attempted to model the shadowed region effect with a combination 
of different ramp models, but nothing we tried was able to cleanly 
separate the effect from the phase variations, and there was always 
some clear structure left behind in the residuals of the fit. Another 
diagnostic that our detrending attempts were unsuccessful was that 
the phase offset as a function of wavelength smoothly varied around 
~10° Ein the unaffected region of the detector, but in the shadowed 
region, the phase offset would abruptly change to -5? W; sucha sharp 
change in a suspect region ofthe detector seems highly unlikely to be 
astrophysical in nature. As a result, we ultimately chose to exclude the 
three spectral bins spanning 10.5-12 jum from our retrieval efforts. 


Determining the number of sinusoid harmonics. To determine the 
complexity of the phase-curve model required to fit the data, we used 
the Eureka! v1 reduction and most of the Eureka! v1 fitting methods 
described below, with the exception of using the dynesty? nested 
sampling algorithm (which computes the Bayesian evidence, z) and 
abatmantransit and eclipse model. Within dynesty, we used 256 live 
points, ‘multi’ bounds, 'rwalk' sampling, and ran until the estimated 
dIn(Z)reached 0.1. We then evaluated first-, second- and fourth-order 
models for the broadband light curve, excluding all third-order sinu- 
soidalterms from the fourth-order model as these terms are not likely 
to be produced by the planet's thermal radiation???', We then com- 
paredthe Bayesian evidences ofthe different models following refs. 
52,53 and found that the second-order model was significantly pre- 
ferred over the first-order model at 160 (AIn(Z) = 128), while the 
second-order model was insignificantly preferred over the 
fourth-order model at 2.20 (AIn(Z) = 1.3). This is also confirmed by 
eye where the first-order model leaves clear phase-variation signa- 
tures in the residuals, while the residuals from the second-order 
model leave no noticeable phase variations behind. Finally, we also 
compared the phase-resolved spectra obtained from different order 
phase-curve models; we found that our spectra significantly changed 
going from a first- to second-order model (altering one or more spec- 
tral points by >10), but the fourth-order model did not significantly 
change the resulting phase-resolved spectra compared with the 
second-order. As a result, the final fits from all reductions used a 
second-order model. Thebroadband light curves obtained from the 
four reductions and the associated phase-curve models are shown in 
Supplementary Fig. 1. 


Eureka! v1 fitting methods. We first sigma-clipped any data points 
that were 40 discrepant from a smoothed version of the data (made 
using a boxcar filter with a width of 20 integrations) to remove any 
obviously errant data points while preserving the astrophysical signals 
like the transit. 

Our astrophysical model consisted of astarry™ transit and eclipse 
model, as well as a second-order sinusoidal phase-variation model. The 
complete astrophysical model had the form 


A(t) = F.(t) + FayE(O (9). (1) 


where tis the time, F.is the received stellar flux (and includes the starry 
transit model), F,;, is the planetary flux at mid-eclipse, E(t) is the starry 
eclipse model (neglecting eclipse mapping signals for the purposes of 
this paper), and V($) is the phase-variation model of the form 


V($) = 1+ AmpCosl x (cos(@) — 1) + AmpSinl x sin(@) (2) 


+AmpCos2 x (cos(2@) — 1) + AmpSin2 x sin(29), 


where @ is the orbital phase in radians with respect to eclipse, and 
AmpCos1, AmpSin1, AmpCos2 and AmpSin2 are all fitted coefficients. 
The second-order phase-variation terms allow for thermal variations 
across the face of the planet that are more gradual or steep than a simple 
first-order sinusoid would allow. We numerically computed dayside, 
morning, nightside and evening spectra using the above U/($) func- 
tion at $ = 0, 1/2, mand 31/2, respectively. To allow the starry eclipse 
function to account for light travel time, we used a stellar radius (R.) of 
0.667 R, (ref. 55) to convert the fitted a/R. (the scaled semi-major axis) 
to physical units. For our transit model, we used a reparameterized 
version of the quadratic limb-darkening model” with coefficients u, 
and u, uniformly constrained between 0 and 1, and used a minimally 
informative prior on the planet-to-star radius ratio (R,/R.). 

Our systematics model consisted of a single exponential ramp in 
time to account for the idle-recovery drift documented for MIRI/LRS 
time-series observations”, a linear trend in time, and a linear trend 
with the spatial position and PSF width. The full systematics model 
can be written as 


SCE, y, Sy) = L(t) x RG) x YCy) x SY(s,), (3) 


Thelinear trend in time is modelled as 


L(tj) = Co + Gb, (4) 


where ¢ is the time with respect to the mid-point of the observations 
and where c, and c; are coefficients. The exponential ramp is modelled 
as 


R(t) 21-4 roe" (5) 


where t, is the time with respect to the first integration and where ry and 
r,are coefficients. The linear trends as a function of spatial position, y, 
are PSF widths, are modelled as 


Y(y)=1+fy (6) 


and 


SY(sy) = 1+ 8Sy, (7) 


where fand gare coefficients. The linear trends as a function of spatial 
position and PSF width are with respect to the mean-subtracted spatial 
position and PSF width. Finally, we also fitted a multiplier (scatter mut) 
to the estimated Poisson noise level for each integration to allow us 
to account for any noise above the photon limit as well as an incorrect 
value for the gain applied in stage 3. 

With an initial fit to the broadband light curve (5-10.5 jum), we 
assumed a zero eccentricity and placed a Gaussian prior on the planet's 
orbital parameters (period, P; linear ephemeris, t;; inclination, i; and 
scaled semi-major axis, a/R.) based on previously published values for 
the planet”. For the fits to the spectroscopic phase curves, wethen fixed 
these orbital parameters to the estimated best fit from the broadband 
light curve fit to avoid variations in these wavelength-independent 
values causing spurious features in the final spectra. We fitted the 
observations using the No U-Turns Sampler (NUTS) from PyMC3” 
with 3 chains, each taking 6,000 tuning steps and 6,000 production 
draws with atarget acceptance rate of 0.95. We used the Gelman-Rubin 
statistic? to ensure the chains had converged. We then used the 16th, 
50th and 84th percentiles from the PyMC3 samples to estimate the 
best-fit values and their uncertainties. 
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Eureka! v2 fitting methods. For the second fit made with Eureka!, 
we proceeded very similarly to the Eureka! v1 fit. We clipped the light 
curves using a boxcar filter of 20 integrations wide with a maximum 
of 20 iterations and a rejection threshold of 4cto reject these outliers. 
We also modelled the phase curve using a second-order sinusoidal 
function, but we modelled the transit and eclipse using batman” 
instead of starry. Like in the Eureka! v1 fit, we modelled instrumental 
systematics with a linear polynomial model in time (equation (4)), an 
exponential ramp (equation (5)), a first-order polynomial in y posi- 
tion (equation (6)) and a first-order polynomial in PSF width in the s, 
direction (equation (7)). 

We fitted the data using the emcee sampler” instead of NUTS, with 
500 walkers and 1,500 steps. The jump parameters that we used were 
the same as in the Eureka! v1 fit: R,/R., Faay, Uy, uj, AMpCos1, AmpSin1, 
AmpCos2, AmpSin2, Co, C, fo, ff, g and scatter nur (multiplier to the 
estimated Poisson noise level for each integration like in the Eureka! 
v1 fit). We used uniform priors on u and u, from O to1, uniform priors 
onAmpCosl, AmpSin1, AmpCos2, AmpSin2 from -1.5 to 1.5 and broad 
normal priors and all other jump parameters. Convergence, mean 
values and uncertainties were computed in the same way as for the 
Eureka! v1 fit. 


TEATRO fitting methods. To measurethe planet's emission as a func- 
tionoflongitude, we modelled the light curves with a phase-variation 
model,aneclipse model, a transit model and an instrument systematics 
model. The phase-curve model, W(t), consists of two sinusoids: one at 
the planet’s orbital period, P, and one at P/2 to account for second-order 
variations. The eclipse model, F(t), and transit model, T(t), are com- 
puted with the exoplanet“ package that uses the starry package". We 
save the eclipse depth, ô., and normalize E(t) to a value of O during the 
eclipse and 1 out of the eclipse, which we then call E(t). We used pub- 
lished transit ephemerides”, a null eccentricity and published stellar 
parameters”. The planet-to-star radius ratio, R,/R., impact parameter, 
b, and mid-transit time, to, are obtained froma fit to the broadband light 
curve. The systematics model, S(t), is composed ofa linear function to 
account for a downwards trend and an exponential function to account 
for the initial ramp. The full model is expressed as: 


F(t) = (PÀ) — D (te) + Be) x ENO + TE) + S(O) (8) 
W(t) = ay sin2nt/P — by) + cy cos(4nt/P — dy) (9) 


S(t) = as e-5st + cst - ds (10) 


where Y(t.) is the value of V at the mid-eclipse time, t.. 

We fit our model to the data using a Markov chain Monte Carlo 
(MCMC) procedure based on the PyMC3 package” and gradient-based 
inference methods as implemented in the exoplanet package. We set 
normal priors for to, R,/R., the stellar density (p-), ay, by, Cs and d; with 
mean values obtained from an initial nonlinear least-squares fit, anor- 
mal prior for a, with a zero mean, uniform priors for the surface bright- 
ness ratio between the planet's dayside and the star (s), b, c, and dy, 
uninformative priors for the quadratic limb-darkening parameters”, 
and allowed for wide search ranges. We ran two MCMC chains with 
5,000 tuning steps and 100,000 posterior samples. Convergence was 
obtained for all parameters (except in one case where a, was negligible 
and b, was unconstrained). We merged the posterior distributions of 
both chains and used their median and standard deviation to infer final 
values and uncertainties for the parameters. We also verified that the 
values obtained from each chain were consistent. 

For the spectroscopic light-curve fits, we fixed all physical param- 
eters to those obtained from the broadband light-curve fit except the 
surface brightness ratio, s, that sets the eclipse depth, we masked the 
transit part of the light curve, and used a similar procedure. After the 


fits, we calculated the eclipse depth, 6,, as s x (R,/R.)’, and computed 
W(t) for the final parameters, W;(t). The planetary flux is W(t) — U(t.) + ôe- 
We computed the uncertainty on the eclipse depth in two different 
ways: from the standard deviation of the posterior distribution of 
Sx (RR), and from the standard deviation of the in-eclipse points 
divided by yN., where N, is the number of in-eclipse points, and took 
the maximum of the two. To estimate the uncertainty on the planet's 
flux, we computed the Ic interval of Y(t) based on the posterior distri- 
butions ofits parameters, computed the lo uncertainty of d,, and added 
them in quadrature to the uncertainty on the eclipse depth to obtain 
more conservative uncertainties. 

The spectra presented in this paper and used in the combined 
spectra are based on system parameters that were derived froma 
broadband light curve obtained in the 5-12 um range, a transit fit in 
which the stellar mass and radius were fixed, a simpler additive model 
in which the phase curve was not turned off during the eclipse, and an 
MCMC run that consisted in two chains of 10,000 tuning steps and 
10,000 posterior draws. Updated spectra based on system parameters 
derived from the broadband light curve obtained in the 5-10.5 uum 
range, a transit fit that has the stellar density as a free parameter, the 
light-curve model shown in equation (8), and two MCMC chains of 
5,000 tuning steps and 100,000 posterior draws are consistent within 
loat every point with those shown here. As we average four reductions 
and inflate the uncertainties during the retrievals, the impact of these 
updates on our results are expected to be marginal. 


SPARTA fitting methods. We use emcee” to fit a broadband light curve 
with the transit time, eclipse time, eclipse depth, four phase-curve 
parameters (C, and D, for the first-order, and C, and D, for the 
second-order sinusoids), transit depth, a/R., b, flux normalization, 
error-inflation factor, instrumental ramp amplitude (A) and 1/timescale 
(1), linear slope in time (m) with respect to the mean ofthe integration 
times (0, and linear slope with trace y position (c,) as free parameters. 
The best-fit transit and eclipse times, a/R.and bare fixed for the spec- 
troscopic light curves. 

For the spectroscopic fits, we then use emcee to fit the free param- 
eters: the eclipse depth, four phase-curve parameters, error-inflation 
factor, flux normalization, instrumental ramp amplitude and 1/time- 
scale, linear slope with time, and linear slope with trace y position. All 
parameters are given uniform priors. 1/timescale is given a prior of 
5-100 d”, but the other priors are unconstraining. In summary, the 
instrumental model is: 


S-F, (1 + Aexp(—t/T) + Gy + m(t — D). (11) 


while the planetary flux model is: 


Fp = E + C(cos(wt) — 1) + D, sin(wt) + C,(cos(2@t) — 1) + D; sin(2wt), 
(12) 


where Eis the eclipse depth and w = 2r/Pisthe planet’s orbital angular 
frequency. Note that the phase variations were set to be zero during 
eclipse. 


Combining independent spectra. Comparing the phase-resolved 
spectra from each reduction (Supplementary Fig. 2), we see that for 
wavelengths below 10.5 um, the spectra are typically consistent, while 
larger differences arise in the >10.5 um region affected by the shadowed 
region effect. For our final, fiducial spectrum, we decided to use the 
median spectrum and inflated our uncertainties to account for disa- 
greements between different reductions. The median phase-resolved 
spectra were computed by taking the median F,/F. per wavelength. 
The uncertainties were computed by taking the median uncertainty 
per wavelength, and then adding in quadrature the RMS between the 
individual reductions and the median spectrum; this minimally affects 
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the uncertainties where there is minimal disagreement and appreci- 
ably increasesthe uncertainties wherethe larger disagreements arise. 

Each reduction also computed a transmission spectrum (Sup- 
plementary Fig. 2), which appears quite flat (within uncertainties) with 
minimal differences between reductions. WASP-43b is not an excellent 
targetfor transmission spectroscopy, however, and thesetransmission 
spectra are not expected to be overly constraining. 


Atmospheric forward models 

GCMS were used to simulate atmospheric conditions, from which syn- 
thetic phase curves and emission spectra were forward modelled and 
compared with the observations. The GCMs used in this study are listed 
in Supplementary Table 1, and details of each simulation are provided 
in Extended Data Table 1 and the following sections. 


Generic PCM. The Generic Planetary Climate Model (Generic PCM) 
is a three-dimensional global climate model designed for modelling 
the atmosphere of exoplanets and for palaeoclimatic studies. The 
model has been used for the study of planetary atmospheres of the 
Solar System^* 55, terrestrial exoplanets”, mini-Neptunes® and hot 
Jupiters®. The dynamical core solves the primitive equations of meteor- 
ology onaArakawa C grid. The horizontal resolution is 64 x 48 (thatis, 
5.625 x 3.75°) with 40 vertical levels, equally spaced in logarithmic scale 
between 10 Pa and 800 bars. Along with the various parameterizations 
of physical processes described in refs. 64-68, the Generic PCM treats 
clouds as radiatively active tracers of fixed radii. 

The model is initialized using temperature profiles from the radia- 
tive-convective one-dimensional model Exo-REM”. The radiative data 
are computed offline using the out-of-equilibrium chemical profiles 
of the Exo-REM run. We use 27 frequency bins in the stellar channel 
(0.261-10.4 um) and 26 in the planetary channel (0.625-324 um), all 
bins including 16 k-coefficients. We start the model from a rest state (no 
winds), with a horizontally homogeneous temperature profile. Models 
are integrated for 2,000 days, which is long enough to complete the 
spin-up phase of the simulation above the photosphere. We do not 
include Rayleigh drag in our models. The simulations are performed 
including clouds of Mg,SiO,, with varying cloud radii (0.1, 0.5, 1, 3 and 
5 um). We also computed cloudless and Mg,SiO, models with a 10x 
solar metallicity and the same radii for the cloud particles. Regardless 
of the composition and size of the clouds, our model clearly indicates 
that there is no cloud formation on the dayside. Asymmetric limbs are a 
natural result of our model, with the eastern terminator being warmer 
while the western limb is cloudier and cooler. Spectral phase curves 
were produced using the Pytmosph3R code”. 


SPARC/MITgcm with radiative transfer post-processing by 
gCMCRT. SPARC/MITgcm couples a state-of-the-art non-grey, 
radiative-transfer code with the MITgcm™. The MITgcm solves the prim- 
itive equations of dynamical meteorology on a cubed-sphere grid”. 
It is coupled to the non-grey radiative-transfer scheme based on the 
plane-parallel radiative-transfer code of ref. 73. The stellar irradiation 
incident on WASP-43b is computed with a PHOENIX model/'*75, We use 
previously published opacities”, including more recent updates/?^, 
and the molecular abundances are calculated assuming local chemical 
equilibrium®®. In the GCM simulations, the radiative-transfer calcula- 
tions are performed on 11 frequency bins ranging from 0.26 um to 
300 um, with 8 k-coefficients per bin statistically representing the 
complex line-by-line opacities’. The strong visible absorbers TiO 
and VO are excluded in our k tables similar to our previous GCMs of 
WASP-43b?? that best match the observed dayside emission spectrum 
and photometry. 

Clouds in the GCM are modelled as tracers that are advected by 
the flow? and can settle under gravity. Their formation and evapora- 
tion are subjected to chemical equilibrium predictions, that is, the 
condensation curves of various minerals described in ref. 80. The 


conversion between the condensable ‘vapour’ and clouds is treated as 
asimple linear relaxation over a short relaxation timescale of 100 s. The 
scattering and absorption of the spatial- and time-dependent clouds 
areincluded in both the thermal and visible wavelengths of the radia- 
tive transfer. A similar dynamics-cloud-radiative coupling has been 
developed in our previous GCMs with simplified radiative transfer and 
has been used to study the atmospheric dynamics of brown dwarfs”? 
and ultrahot Jupiters*?. Clouds are assumed to follow a log-normal 
size distribution?*, which is described by the reference radius ry and a 
N [Inr/ro) 
V2nor exp ( 202 
the number density per radius bin of r and v is the total number 
density. o and r, are free parameters and the local 3 is obtained from 
the local mass mixing ratio of clouds. The size distribution is held fixed 
throughout the model and is the same for all types of cloud. 

Our GCMs do not explicitly impose a uniform radiative heat flux 
at the bottom boundary but rather relax the temperature of the low- 
est model layer (that is, the highest pressure layer) to a certain value 
over a short timescale of 100 s. This assumes that the deep GCM layer 
reaches the convective zone and the temperature there is set by the 
interior convection that ties to the interior structure of the planet. This 
lowest-layer temperature is in principle informed by internal structure 
models of WASP-43b, which are run by MESA hotJupiter evolution 
modules" to match the present radius of WASP-43b. In most models, 
this lowest-layer temperature is about 2,509 K at about 510 bars. The 
horizontal resolution of our GCMsis typically C48, equivalent to about 
1.88? per grid cell. The vertical domain is from2 x 10 * bar at the top to 
700 barsatthebottom andis discretized to 53 vertical layers. We typi- 
cally run the simulation for over 1,200 days and average all physical 
quantities over the last 100 days of the simulations. 

All our GCMs assume a solar composition. We performed a base- 
line cloudless model and one case with only MnS and Na;S clouds with 
ro = 3 um, and then a few cases with MnS, Na,S and MgSiO, clouds with 
ry71,1.5,2and3 um. The ois held fixed at 0.5 in all our cloudy GCMs. 

We post-process our GCM simulations with the state-of-the-art 
gCMCRT code, which is a publicly available hybrid Monte Carlo radi- 
ative transfer (MCRT) and ray-tracing radiative-transfer code. The 
model is described in detail in ref. 85 and has been applied to a range 
ofexoplanet atmospheres®*. gCMCRT can natively compute albedo, 
transmission and emission spectra at both low and high spectral reso- 
lution. gCMCRT uses custom k tables, which take cross-section data 
from both HELIOS-K* and EXOPLINES**. Here, we apply gCMCRT to 
compute low-resolution emission spectra and phase curves at R « 300 
from our GCM simulations. We use the three-dimensional temperature 
and condensate cloud tracer mixing ratio from the time-averaged 
end-state of each case. We assume the same cloud particle size distri- 
bution as our GCMs. 


non-dimensional deviation o: n(r) 


) wheren(r)is 


expeRT/MITgcm. The GCM expeRT/MITgcm uses the same dynami- 
cal core as SPARC/MITgcm and solves the hydrostatic primitive 
equations on a C32 cubed-sphere grid”. It resolves the atmosphere 
above 100 bar on 41 log-spaced cells between 1 x 10? bar and 100 bar. 
Below100 bar, 6 linearly spaced grid cells between 100 bar and 700 bar 
are added. The model expeRT/MITgcm thus resolves deep dynamics 
in non-inflated hot Jupiters like WASP-43b'^9?, 

TheGCM is coupledto a non-grey radiative-transfer scheme based 
on petitRADTRANS”. Fluxes are recalculated every fourth dynamical 
time step. Stellar irradiation is described by the spectral fluxes from the 
PHOENIX model atmosphere suite” *. The GCM operates ona precalcu- 
lated grid of correlated k-binned opacities. Opacities from the ExoMol 
database” are precalculated offline ona grid of 1,000 logarithmically 
spaced temperature points between 100 K and 4,000 K for every ver- 
tical layer. We further include the same species as shown in ref. 89 
except TiO and VO to avoid the formation of a temperature inversion 
in the upper atmosphere. These are: H,O (ref. 92), CH, (ref. 93), CO, 
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(ref. 94), NH, (ref. 95), CO (ref. 96), H,S (ref. 97), HCN (ref. 98), PH, 
(ref. 99), FeH (ref. 100), Na (refs. 74,101) and K (refs. 74,101). For Ray- 
leigh scattering, the opacities are H, (ref. 102) and He (ref. 103), and we 
add the following collision-induced absorption (CIA) opacities: H.-H, 
(ref. 104) and H;-He (ref. 104). We use for radiative-transfer calcula- 
tions in the GCM the same wavelength resolution as SPARC/MITgcm 
(S1), but incorporate 16 instead of 8 k-coefficients. Two cloud-free 
WASP-43b GCM simulations were performed, one with solar and one 
with 10x solar element abundances. Each simulation ran for 1,500 days 
to ensure that the deep wind jet has fully developed. The GCM results 
used inthis paper were time averaged over the last 100 simulation days. 

Spectra and phase curves are produced from our GCM results in 
post-processing with petitRADTRANS” and prt_phasecurve® using a 
spectral resolution of R = 100 for both the phase curve and the spectra. 


RM-GCM. Originally adapted from the GCM of ref. 105 by refs. 106-108, 
the RM-GCM model has been applied to numerous investigations of hot 
Jupiters and mini-Neptunes?*'^?"!, The GCM's dynamical core solves 
the primitive equations of meteorology using a spectral representa- 
tion of the domain, and it is coupled to a two-stream, double-grey 
radiative-transfer scheme based on ref. 112. Recent updates have added 
aerosol scattering? with radiative feedback?^?*. Following ref. 8, aero- 
solsare representative of condensate clouds and are treated as purely 
temperature-dependent sources of opacity, with constant mixing 
ratios set by the assumed solar elemental abundances. The optical 
thicknesses of the clouds are determined by converting the relative 
molecular abundances (or partial pressures) of each species into par- 
ticles with prescribed densities and radii?. The model includes up to 13 
different cloud species of various condensation temperatures, abun- 
dances and scattering properties. Places where clouds overlap have 
mixed properties, weighted by the optical thickness of each species. 

Simulations from this GCM included a clear atmosphere and two 
sets of cloudy simulations. Following ref. 8, one set of cases included 13 
different species: KCI, ZnS, Na,S, MnS, Cr,0;, SiO}, Mg,SiO,, VO, Ni, Fe, 
Ca,SiO,, CaTiO, and Al,O,; the other set omitted ZnS, Na,S, MnS, Fe and 
Ni, based on considerations of nucleation efficiency’’. For both cloud 
composition scenarios, the models explored the observational conse- 
quences of variations in the cloud deck’s vertical thickness through a 
series of simulations with clouds tops truncated over a range of heights 
at 5-layer intervals (roughly a scale height), ranging from 5to 45 layers 
ofthe 50-layer model. This effectively mimics a range of vertical mix- 
ingstrengths. From the complete set published in ref. 26, we selected 
a subset, with clouds of maximum vertical extents between two and 
nine scale heights from each ofthe two cloud composition scenarios. 

Simulations were initialized with clear skies, no winds and no hori- 
zontal temperature gradients. We ran the simulations for over 3,500 
planetary orbits, assuming tidal synchronization. Resulting tempera- 
ture, wind and cloud fields of the GCM were then post-processed''*!^ 
to yield corresponding emission phase curves. 


THOR. THOR*'" is an open-source GCM developed to study the 
atmospheres and climates of exoplanets, free from Earth- or Solar 
System-centric tunings. The core that solves the fluid flow equations, 
the dynamical core, solves the non-hydrostatic compressible Euler 
equations on an icosahedral grid"°"*. THOR has been validated and 
used to simulate the atmosphere of Earth”, Solar System planets”?! 
and exoplanets!6177722 

For this work, THOR used the same configuration as with previ- 
ously published simulations to study the atmospheric temperature 
structure, cloud cover and chemistry of WASP-43b^^*!?, Two simula- 
tions were conducted, one with a clear atmosphere and another with a 
cloud structure on the nightside ofthe planet. To represent the radia- 
tive processes, THOR uses a simple two-band formulation calibrated 
to reproduce the results from more complex non-grey models on 
WASP-43b*"*, A simple cloud distribution on the nightside of the planet 


and optical cloud properties are parameterized‘ and adapted to repro- 
duce previous HST” and Spitzer*” observations. These simulations on 
WASP-43b with THOR have also been used to test the performance of 
future Ariel phase-curve observations”. 

Both simulations, with clear and cloudy atmospheres, started 
with isothermal atmospheres (1,440 K, equilibrium temperature) and 
integrated for roughly 9,400 planetary orbits (assuming atidally locked 
configuration) until a statistically steady state of the deep atmosphere 
thermal structure was reached. The long integration avoids biasing the 
results towards the set initial conditions’”°. 

The multiwavelength spectra are obtained from post-processing 
the three-dimensional simulations with a multiwavelength 
radiative-transfer model’. The disk-averaged planet spectrum is 
calculated at each orbital phase by projecting the outgoing inten- 
sity for each geographical location of the observed hemisphere. The 
spectra include cross-sections of the main absorbers in the infrared, 
drawn from the ExoMOL (H,O (ref. 92), CH, (ref. 127), NH, (ref. 128), 
HCN (ref. 129) and H,S (ref. 97)), HITEMP?? (CO, and CO) and HITRANP! 
(C;H;) databases. The Na and K resonance lines’” are also added, as 
were H;-H, and H,-He CIA"*, The atmospheric bulk composition 
was assumed to have solar abundance (consistent with HST/WFC3 
spectrum observations), and each chemical species concentration 
was calculated with the FastChem model”. The PHOENIX models” 
were used for the WASP-43 star spectrum. 


Atmospheric retrieval models 

We perform atmospheric retrievals on the phase-resolved emission 
spectra using six different retrieval frameworks, each described in 
turn below. The chemical constraints from these analyses are sum- 
marized in Extended Data Tables 2 and3, and the spectral fits obtained 
are shown in Extended Data Fig. 3. Across the six retrieval analyses, we 
use anerror-inflation parameter to account for the effects of unknown 
data and/or model uncertainties. This free parameter is wavelength 
independent and multiplies the 1c error bars in the calculation of the 
likelihood function in the Bayesian sampling algorithm. 


HyDRA retrieval framework. The HyDRA atmospheric retrieval 
framework" consists of a parametric atmospheric forward model 
coupled to PyMultiNest "5, a nested sampling Bayesian parameter 
estimation algorithm". HyDRA has been applied to hydrogen-rich 
atmospheres”, and further adapted for secondary atmospheres”? 
and high-resolution spectroscopy in both one and two dimensions", 
The input parameters for the atmospheric forward model include 
constant-with-depth abundances for each of the chemical species 
considered, six temperature profile parameters corresponding tothe 
temperature profile model of ref. 143, and aconstant-with-wavelength 
multiplicative error-inflation parameter to account for model uncer- 
tainties. We additionally include a dilution parameter, Aus, for the 
dayside, morning and evening hemispheres, which multiplies the emis- 
sion spectrum by a constant factor «1 and accounts for temperature 
inhomogeneities in each hemisphere'**. 

We consider opacity contributions from the chemical species 
that are expected to be present in hot Jupiter atmospheres and that 
have opacity in the MIRI LRS wavelength range: H,O (ref. 130), CH, 
(refs. 127,145), NH; (ref. 128), HCN (refs. 98,129), CO (ref. 130), CO, 
(ref. 130), C,H, (refs. 131,146), SO, (ref. 147), H,S (refs. 97,148) and 
CIA due to H,-H, and H,-He (ref. 104). The line-by-line absorption 
cross-sections for these species are calculated following the methods 
described in ref. 134, using data from each of the references listed. We 
further explore retrievals with a simple silicate cloud model, which 
includes the modal particle size, cloud particle abundance, cloud 
base pressure and a pressure exponent for the drop-off of cloud 
particle number density with decreasing pressure. The opacity struc- 
ture of the cloud is calculated using the absorption cross-sections 
of ref. 149. 
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Given the input chemical abundances, temperature profile and 
cloud parameters, the forward model calculates line-by-line radiative 
transfer to produce the thermal emission spectrum at a resolution 
of R «15,000. The spectrum is then convolved to a resolution of 100, 
binned tothe same wavelength bins as the observations and compared 
with the observed spectrum to calculate the likelihood of the model 
instance. The nested sampling algorithm explores the parameter space 
using 2,000 live points, and further calculates the Bayesian evidence of 
the retrieval model, which can be used to compare different models”. 
In particular, we calculate the detection significance of a particular 
chemical species by comparing retrievals that include/exclude that 
species, fixing the value of the error-inflation parameter to be the 
median retrieved value found with the full retrieval model. 

Across the four phases, the only chemical species detected with 
statistical significance (230) is H,O. The retrieved H,O abundances 
are in the range ~30-10* ppm (1c uncertainties), with detection sig- 
nificances varying between -30 and -4o (Extended Data Table 2). We 
do not detect CH, at any phase, and place an upper limit of 16 ppm on 
the nightside CH, abundance, potentially indicating disequilibrium 
chemistry processes as described in the main text. We do not detect 
NH, at any phase either, consistent with the very low NH, abundances 
predicted by both chemical equilibrium and disequilibrium models”. 
The retrievals do not favour cloudy models over clear models with 
statistical significance, with extremely weak preferences of «1o at all 
phases. In addition, the posterior probability distributions for the 
cloud parameters are unconstrained. Extended Data Fig. 5 shows the 
pressure ranges ofthe atmospheric model probed by the observations. 


PyratBay retrieval framework. PyratBay is an open-source framework 
that enables atmospheric modelling, spectral synthesis, and atmos- 
pheric retrievals of exoplanet observations ?. The atmospheric model 
consists of parametrictemperature, composition and altitude profiles 
as a function of pressure, for which emission and transmission spec- 
tra can be generated. The radiative-transfer module considers opac- 
ity from alkali lines?', Rayleigh scattering’”"’, Exomol and HITEMP 
molecular line lists? ^* pre-processed with the REPACK package"? 
to extract the dominant line transitions, CIAP? and cloud opacities. 
The PyratBay retrieval framework has the ability to stagger model 
complexity and explore a hierarchy of different model assumptions. 
Temperature models range from an isothermal profile to physically 
motivated parameterized models", Composition profiles range 
from the simpler constant-with-altitude ‘free abundance’ to the more 
complex ‘chemically consistent’ retrievals, the latter accomplished 
via the TEA code’; while cloud condensate prescriptions range from 
the classic ‘power law + grey’ to a ‘single-particle-size’ haze profile, a 
partial-coverage factor ‘patchy clouds'??, and the complex parameter- 
ized Mie-scattering thermal stability cloud (TSC) model (J.B. et al., man- 
uscript in preparation). The TSC cloud prescription, initially inspired 
by refs. 84,160, has additional flexibility in the location of the cloud 
base and was further improved for this analysis (see below). The for- 
mulation utilizes a parameterized cloud shape, effective particle size 
and gas number density below the cloud deck, while the atmospheric 
mixing and settling are wrapped up inside the cloud extent and the 
condensate mole fraction as free parameters. This cloud model was 
applied to WASP-43b JWST/MIRI phase-curve simulations”, gener- 
ated during the JWST preparatory phase, in anticipation of the actual 
WASP-43b JWST/MIRI observations. We showed that the TSC model 
has the ability to distinguish between MgSiO, and MnS clouds on the 
nightside ofthe planet. 

For this analysis, we conducted a detailed investigation using 
various model assumptions. We started by exploring simple tem- 
perature prescriptions and gradually moved towards more complex 
ones. Initially, we considered opacity contributions from all chemical 
species expected to be observed in the MIRI wavelength range (H,O, 
CH, NH,, HCN, CO, CO,, C;H,, SO,, H5S), but eventually focused on 


only those that are fit by the data. We also implemented the dilution 
parameter" and an error-inflation factor, which account for some 
additional model and data uncertainties. The constraints on H,O 
(together with the detection significance?) and the upper limit for 
CH, for all phases are given in Extended Data Table 2. The abundances 
of these species across all phases were largely model independent. 
However, the tentative constraints on NH,, which we saw in multiple 
phases, werestrongly model dependent, and were completely erased 
with the inclusion of the dilution parameter and the error inflation, 
thus we do not report them here. WASP-43b emission spectra were 
computed at a resolution of R « 15,000 utilizing opacity sampling of 
high-resolution pre-computed cross-sections (R ~ 10°) of considered 
species. Furthermore, we thoroughly examined the possibility of 
detecting clouds in each of the four-quadrant phases, with a special 
emphasis on the nightside of the planet. To do this, we employed 
the TSC model, as in our previous analysis”, and explored a range of 
cloud species, MgSiO,, MnS, ZnS and KCI, that would condense under 
the temperature regimes expected for WASP-43b!? (Extended Data 
Fig. 6, left). We also introduced the effective standard deviation of 
the log-normal distribution? as a free parameter (Oog), allowing for 
even more flexibility in our cloud model (Extended Data Fig. 6, right, 
last subpanel). To thoroughly explore the parameter space, we used 
two Bayesian samplers, the differential-evolution MCMC algorithm! 
implemented following ref. 164, and the nested sampling algorithm, 
implemented through PyMultiNest’>””®, utilizing 15 million models 
and 2,000 live points, respectively. Our investigation did not provide 
constraints on any of the cloud parameters for any of the explored 
cloud condensates at any of the planetary phases, indicating the 
absence of detectable spectral features from clouds in the observa- 
tions (Extended Data Fig. 6, right). 


NEMESIS retrieval framework. NEMESIS 9*5 is a free retrieval frame- 
work that uses a fast correlated-K'*' forward model, combined with 
either an optimal estimation or nested sampling retrieval algorithm. 
It has been used to perform retrievals on spectra of numerous plan- 
etary targets, both inside and outside the Solar System’. In this 
work, we use the PyMultiNest sampler? with 500 live points. The 
retrieval model presented includes four spectrally active gases, H,O 
(ref. 92), CO (ref. 96), CH, (ref. 93) and NH, (ref. 95), with k tables cal- 
culated as in ref. 91; we did not include CO, or H.S after initial tests 
indicated these were not required to fit the spectrum. All gases are 
assumed to be well mixed in altitude. CIA from H, and He is taken from 
refs. 156,170. The spectrum is calculated at the resolution of the obser- 
vation, using optimized channel integrated k tables generated from 
original k tables with a resolving power R = 1,000. The temperature 
profile is modelled as athree-parameter Guillot profile, after ref. 157, 
with free parameters x, y and £ (a is fixed to be zero). We include a 
well-mixed, spectrally grey cloud with a scalable total optical depth 
with a cloud top at 12.5 mbar. The other retrieved parameters are 
a hotspot dilution factor for phases 0.25, 0.5 and 0.75, following 
ref. 144, and an error-inflation term. 

Tocalculatethe detection significance for H,O, we run the retrieval 
with and without HO, with all other aspects of the run identical. We 
then take the difference of the PyMultiNest global log-evidence val- 
ues for the two scenarios, and convert from log(Bayesian evidence) 
to sigma following ref. 52. The 99% upper limit for CH, is calculated 
from the equally weighted posterior distribution. We also attempt to 
retrieve CO and NH, abundances. CO is generally poorly constrained, 
and NH, is unconstrained for phases O and 0.75; for log(NH;), we recover 
a99% upper limit of 22.2 at phase 0.25 and -3.9 at phase 0.5. The cloud 
opacity is also generally unconstrained, with the total optical depth 
able to span several orders of magnitude. We stress that this model 
is very crude as it has only one variable cloud parameter, and further 
exploration of suitable cloud models for mid-infrared phase curves is 
warranted in future work. 
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SCARLET retrieval framework. We perform atmospheric retrievals 
onthe four phase-resolved spectra using the SCARLET framework”, 
The planetary disk-integrated thermal emission, F,, is modelled for a 
given set of atomic/molecular abundances, temperature- pressure 
profile and cloud properties. We compare our model spectra with the 
observations by normalizing the thermal emission F, using a PHOE- 
NIX” ” stellar model spectrum with effective temperature Tır = 4,300 K 
and surface gravity log g = 4.50. The model spectra are computed ata 
resolving power of R = 15,625, convolved to the resolving power of MIRI/ 
LRS and then binned to the 11 spectral bins («10.5 um) considered in 
the analysis, assuming the throughput to be uniform over asinglebin. 

The atmospheric analysis is performed considering thermochemi- 
cal equilibrium, where the metallicity [M/H] (w[—3, 3]) and carbon-to- 
oxygen ratio (u[0, 3) are free parameters that dictate the overall atmos- 
pheric composition. We use a free parameterization of the tempera- 
ture-pressure profile’” by fitting for N=4 temperature points 
(u[100, 4400] K) with a constant spacing in log-pressure. The tempera- 
ture-pressure profileis interpolated to the 50 layers (P = 107-10 $ bar) 
considered in the model using a spline function to produce a smooth 
profile. We use a grid of chemical equilibrium abundances produced 
with FastChem2'? to interpolate the abundance of species as a function 
of temperature and pressure for given values of [M/H] and C/O. The 
species considered in the equilibrium chemistry are H, H` (refs. 174,175), 
H, He, H,O (ref. 92), OH (ref. 130), CH, (ref. 127), C,H, (ref. 176), CO 
(ref. 130), CO, (ref. 130), NH, (ref. 95), HCN (ref. 98), PH, (ref. 99), 
TiO (ref. 177) and VO (ref. 178). All opacities for these species are con- 
sidered when computing the thermal emission. We account for poten- 
tial spatial atmospheric inhomogeneities in the planetary disk that are 
observed ata given phase by including an area fraction parameter Aps 
(u[0, 1]), which is meant to represent the possibility of a fraction of the 
disk contributing to most of the observed thermal emission’. This 
parameter is considered for all phases with the exception of the night- 
side, which is expected to be relatively uniform. Finally, we fit for an 
error-inflation parameter k, (u[0.1,10]) to account for potential model 
and data uncertainty, which results in a total of 8 (7 for the nightside) 
free parameters. We consider 8 walkers per free parameter for 
the retrievals which are run for 30,000 steps. The first 18,000 steps 
are discarded when producing the posterior distributions of the 
free parameters. 


PLATON retrieval framework. PLATON'^, Planetary Atmosphere 
Tool for Observer Noobs, is a Bayesian retrieval tool that assumes 
equilibrium chemistry. We adopt the temperature-pressure profile 
parameterization of ref. 180, and use the dynesty nested sampler? to 
retrieve the following free parameters: stellar radius; stellar tempera- 
ture; the log metallicity, log(Z); C/O; 5 temperature-pressure param- 
eters (log(x,,), log(y), log(y2), a, B); and an error multiplier. The stellar 
radius and temperature are given Gaussian priors with means and 
standard deviations set by the measurements in ref. 55: 4,400 + 200 K 
and 0.667 + 0.011 R,, respectively. The combination of the two havea 
similar effect to the dilution parameter of other retrieval codes, which 
multiplies the emission spectrum by a constant. For phase 0.0, we 
obtain a significantly better fit when methane opacity is set to zero 
(thus removing all spectral features from methane). We therefore 
adopt this as the fiducial model, whereas for other phases, we do not 
zero out any opacities. 

For all retrievals, we use nested sampling with 1,000 live points. 
The opacities (computed at R = 10,000) and the line lists used to com- 
pute them are listed in ref. 179. We include all 31 species in retrieval, 
notably including H,O, CO, CO;, CH, (except on the nightside), H;S 
and NH}. 


ARCiS retrieval framework. ARCiS (Artful modelling code for exo- 
planet science) is an atmospheric modelling and Bayesian retrieval 
code?" ? that utilizes the MULTINEST'? Monte Carlo nested sampling 


algorithm. The code was used in previous retrievals ofthe atmosphere 
of WASP-43b in transmission’, using the observations of ref. 184, and in 
phase-resolved emission ^, using the observations of refs. 21,22,25,186. 
Reference /? found some evidence that AIO improves the fit of the 
transmission spectra of WASP-43b in the 1.1-1.6 um region. We there- 
foreincludein our modelsfor this workthe following set of molecules 
inour free molecular retrievals: H,O (ref. 92), CO (ref. 96), CO, (ref. 94), 
NH, (ref. 95), CH, (ref. 93) and AIO (ref. 187). The molecular linelists are 
from the ExoMol^*5* or HITEMP?? databases as specified, and k tables 
from the ExoMolOP opacity database”. CIA for H, and He are taken 
from refs. 156,170. We explore the inclusion of a variety of additional 
molecules that have available line list data with spectral features in the 
region of our observations, including HCN (ref. 98), SiO (ref. 189) and 
N,O (ref. 130). We use the Bayes factor, which is the difference between 
the nested sampling global log-evidence (log E) between two models, 
to assess whether the inclusion ofa particular parameter is statistically 
significant. For this, we run a retrieval with the base set of species only 
and another with the base set plus the molecule being assessed. The 
difference in log E between the two models is converted to a signifi- 
canceinterms of ousing the metric of ref. 52. We explore the inclusion 
of asimple grey, patchy cloud model, which parameterizes cloud top 
pressure and degree of cloud coverage (from O0 for completely clear 
tolfor completely covered). We use 1,000 live points and a sampling 
efficiency of 0.3 in MULTINEST for all retrievals. 

We run retrievals both including and not including a retrieved 
error-inflation parameter. The error-inflation parameter is imple- 
mented as per ref. 190 to account for underestimated uncertainties 
and/or unknown missing forward model parameters. All phases apart 
from 0.0 retrieved a parameter that increases the observational error 
barsby twotothreetimestheiroriginal values. The pressure-tempera- 
ture profile parameterization of ref. 191 is used in all cases. We find 
evidence for the inclusion of H,O for all four phases, although this 
evidence goes from strong to weak when error inflation is included for 
the morning phase (0.75). We find no strong evidence for CH, at any 
phase, with 95% confidence upper limits on the log of the volume mix- 
ing ratio (VMR) of -4.9, -2.9, -3.2 and -2.2 for phases 0.0, 0.25, 0.5 and 
0.75, respectively. We find some model-dependent hints of moderate 
evidence (based on the metric of ref. 52) of 4.40 for NH, at phase 0.5 
(constrained to log VMR = —4.5'92), 3.10 for CO at phase 0.5 
(log VMR = —1.703) and 2.60 for CO at phase 0.25 (log VMR = 4.070). 
However, these disappear when the error-inflation parameter is intro- 
duced. Weare not able to constrain any of the cloud parameters for any 
phase, and so do not finda statistical reason to include our simple cloud 
parameterization in the models to better fit the observations. 


Data availability 

The data used in this paper are associated with JUST DD-ERS pro- 
gramme 1366 (principal investigators N.M.B., J.L.B. and K.B.S.; obser- 
vation 11) and are publicly available from the Mikulski Archive for 
Space Telescopes (https://mast.stsci.edu). Additional intermediate 
and final results from this work are archived on Zenodo at https://doi. 
org/10.5281/zenodo.10525170 (ref. 192). 


Code availability 

We used the following codes to process, extract, reduce and analyse 
the data: STScI's JWST Calibration pipeline*$, Eureka! ^, TEATRO, 
SPARTA*, Generic PCM^^ 55, SPARC/MITgcm???? expeRT/GCM2°?, 
RM-GCM??*106-108. THOR*!9172124426493. HyDRA™, PyratBay’”, NEM- 
ESIS'?65 SCARLET", PLATON™’, starry’, exoplanet?', PyMC3", 
emcee’, dynesty?, numpy™, astropy*"?* and matplotlib’”. 
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Extended Data Fig. 1| The underestimation of uncertainties as a function 
of spectral binning for the L168-9b commissioning observations. a, The 
observed L168-9b transmission spectrum with lo error bars for spectrally 
unbinned data (grey circles), 0.15 zm bins (black squares), 0.5 zm bins (large 
red circles), and a 5-12 um broadband bin (horizontal blue shaded region). The 
spectrum for wavelength pairs is not shown to avoid excessive clutter. b, The 
median of the transit depth uncertainties are shown with blue squares, while 
theobserved scatter in the transmission spectrum is shown with orange circles. 


For unbinned data, the transmission spectrum shows about 2.5 x the scatter 
predicted by the fits to the individuallight curves. Binning pairs of wavelengths 
reduces the level of underestimation of the scatter in the transmission spectrum, 
butconsiderable excess noise remains. Coarser binning schemes like the 
constant 0.15 m bins used in the MIRI time-series observation commissioning 
paper” or the 0.5 um bins we use in this work further reduce the level of 
uncertainty underestimation. 
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Extended Data Fig. 2 | A model-independent demonstration of the initial 
changes in flux for the WASP-43b observations. a, The first 120 minutes 

of three of our spectroscopically binned light curves of WASP-43b (with 10 
uncertainties) showing the initial settling behaviour as a function of wavelength. 
Ateal dashed line shows the amplitude of a -0.25% change in flux compared to the 
values around 120 minutes, and a magenta dotted line shows a +0.25% change. 

b, Asummary of the ramp amplitudes, signs, and timescales for each of our 
wavelength bins (with lo uncertainties). The teal and magenta horizontal lines are 
the same as those in panel a to aid in translating between the two figures. At short 
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wavelengths, the flux sharply drops by about 0.5% within the first 30 minutes 
andthen largely settles but does continue to decrease with time. With increasing 
wavelength, the strength of this initial ramp decreases and eventually changes 
sign, becoming an upwards ramp. Within the 'shadowed region' (marked in red), 
thelight curves show a very strong upwards ramp that takes much longer (greater 
than about 60 minutes) to appreciably decay. It is important to note that the 

data in this figure also includes a small amount of astrophysical phase variations 
which should result in a small increase in flux of less than 0.05% per hour. 
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Extended DataFig. 3 | Retrieved spectra from the six retrievals. a, Median opacity. b, c, and d, Same as panel a for the evening terminator, dayside, and 
retrieved nightside spectra for the HyDRA (dark blue line), NEMESIS (dash- morning terminator respectively. e, f, g, and h, Same as panels a, b, c, and d, for 
dotted gold line), and PyratBay (dashed magenta line) and their 1o contour. The the SCARLET (dashed red line), PLATON (blue line), and ARCiS (dash-dotted 
regions of higher water opacity are indicated by the purple shading at the top green line) retrievals. 


of the panel, with the observed rise in flux at 6.3 rm being caused by a drop in 
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Extended Data Fig. 4 | Chemically-consistent atmospheric retrievals. 

Same as Figure 4 but for retrievals assuming thermochemical-equilibrium 
abundances consistent with the pressure-temperature profiles. a, 1o credible 
interval contours ofthe temperature profiles. The black curves show the 
predicted temperature profile from a 2D radiative-transport model46. The 
vertical bars show the range of pressures probed by the observations. b and 

c, probability posterior distributions for H,O and CH, abundances, respectively. 


3000 
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The shaded area for each curve denotes the lo credible interval of each posterior. 
Thegreen and blue bars denote the abundances predicted by equilibrium 

and disequilibrium-chemistry models with solar abundances, respectively, at 
the pressures probed by the observations. Compared to the free-chemistry 
retrievals, the thermochemical-equilibrium retrievals on the nightside spectra 
produced worse fits, this is driven particularly by the higher amount of methane 
expected under equilibrium chemistry. 
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Extended Data Fig. 5 | Retrieval contribution functions. Contribution from the water band around 7-9 utm makes these wavelengths probe lower 
functions integrated over the data point spectral bins, at each phase (a-d), and pressures and hence colder temperatures, whereas the rest of the observing 
for each retrieval framework. These curves show the range of pressures probed window probes higher pressures and higher temperatures. 


by the observation according to the atmospheric models. The enhanced opacity 
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Extended Data Fig. 6 | PyratBay clouds exploration. a, Cloud species that 
condense in the temperature regime expected for the WASP-43b nightside. 
Dashed lines represent vapour pressure curves'? for each species assuming solar 
composition, while the coloured ranges denote the corresponding extent ofthe 
vapour pressure curves assuming 100 x sub- and super-solar atmospheric 
composition. The extent ofthe retrieved nightside contribution functions is 
shown in grey, and the extent of the retrieved temperature uncertainties is shown 
inlight purple. The intersection between the contribution function and 
temperature ranges indicates the pressures at which we could observe cloud 
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condensation and potentially detect their spectral features, if presentin the 
observations. b, Panels display the retrieved posterior density plots for the 
explored cloud parameters of the TSC model (cloud number density, q*; effective 
particle size, re; and the standard deviation of the log-normal distribution, diog) 
for the MnS clouds. The black vertical line denotes the parameter’s median value, 
while the extent of the purple region denotes the Ic uncertainties, both given at 
the top left corner ofthe panel. Similar, fully non-constrained posteriors are 
retrieved for other explored cloud species, MgSiO,, ZnS, and KCl, suggesting the 
lack of observable spectral characteristics from clouds in the observed data. 
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Extended Data Fig. 7 |Acomparison of the retrieved temperature-pressure viewing angle, to produce a one-dimensional profile that is comparable to the 
profiles to the GCM simulations. Each of a-d shows the temperature profile retrieved profile. The GCM simulations are generally warmer on the nightside 
retrieved by HyDRA, compared to the GCM simulations highlighted in Figure3 thanthe retrieved temperatures; cloudy simulations emit from lower pressures 
and listed in Extended Data Table 1. The GCM temperature profiles are calculated and so match the observed lower brightness temperatures better (see the 
at phases 0.0, 0.25, 0.5, and 0.75 by averaging over the visible hemisphere by contribution functions in Extended DataFig. 5). 
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Extended Data Table 1| General Circulation Model simulations used in this investigation 


GCM Cloud Model Metallicity Phase Curve: 
(xsolar) max min offset 
1.* Generic PCM Cloudless 1 5603 2898 31.0 
2. Generic PCM Cloudless 10 5900 2246 16.6 
3.* Generic PCM Mg»SiO,;, 0.1 um 1 5081 998 16.2 
4. Generic PCM Mg»?SiO,, 0.1 um 10 5800 744 9.7 
5. Generic PCM Mg»SiO,, 0.5 um 1 5246 965 18.0 
6. Generic PCM Mg»SiO,;, 0.5 uum 10 5807 697 10.1 
7. Generic PCM Mg2SiOu, 1 um 1 5470 925 16.2 
8. Generic PCM Mg2SiO4, 1 um 10 5816 670 10.8 
9. Generic PCM Mg2SiOu, 3 um 1 5563 990 23.0 
10. Generic PCM Mg2SiOu, 3 um 10 5864 810 14.0 
11. Generic PCM Mg2Si0O,4, 5 um 1 5586 1152 25.9 
12. Generic PCM Mg»?SiO,;, 5 um 10 5880 910 14.8 
13.* SPARC/MITgcm Cloudless 1 5833 2732 29.9 
14. SPARC/MITgcm MnS, Na5S, MgSiOs, 1 um 1 5715 543 11.5 
15. SPARC/MITgcm MnS, Na5S, MgSiOs, 1.5 um 1 5725 533 15.1 
16. SPARC/MITgcm MnS, Na2S, MgSiOs, 2 um 1 5785 638 17.3 
17:* SPARC/MITgcm MnS, Na5S, MgSiOs, 3 um 1 5987 1914 28.8 
18. SPARC/MITgcm MnS, Na5S, 3 um 1 5883 2170 29.2 
19. expeRT/GCM Cloudless 1 5504 2616 34.6 
20.* expeRT/GCM Cloudless 10 5808 2158 19.4 
21.* RM-GCM Cloudless 1 5974 2681 28.1 
22. RM-GCM 13 cloud species, 2 H 1 6048 1828 19.1 
23. RM-GCM 13 cloud species, 3 H 1 5839 1638 16.2 
24. RM-GCM 13 cloud species, 7 H 1 6974 1441 16.9 
25. RM-GCM 8 cloud species, 2 H 1 6007 1919 16.6 
26.* RM-GCM 8 cloud species, 3 H 1 5843 1395 11.2 
27. RM-GCM 8 cloud species, 5 H 1 5109 1040 10.4 
28. RM-GCM 8 cloud species, 7 H 1 6520 1421 18.7 
29. RM-GCM 8 cloud species, 9 H 1 7274 1503 15.5 
30.* THOR Cloudless 1 5248 3471 19.1 
31.* THOR Fixed nightside grey clouds 1 5474 1843 7.6 


For cloud models, the composition of the condensates are provided, along with the prescribed mean radius of the condensed cloud particles. For the RM-GCM cloud models, simulations 
with eight species of clouds include KCl, Cr,O3, SiO;, Mg,SiO,, VO, Ca;SiO,, CaTiO,, and ALO,, while those with 13 species additionally include ZnS, Na;S, MnS, Fe, and Ni; in these cases, 
particle sizes vary with height, and the prescribed maximum vertical extent of the clouds, expressed in pressure scale heights (H), are also provided. Further details of the models can be 
found in the text. Phase curve maxima and minima are expressed as the planet-to-star flux ratio in ppm, while offsets are expressed in degrees of longitude east relative to the sub-stellar 
longitude. Simulations featured in Figure 3 are indicated (*). 
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Extended Data Table 2 | Inferences on H,O and CH, abundances and detection significances 


log(Xu5o) (detection significance) 
Phase: 0.0 0.25 0.5 0.75 


HyDRA | -3.5414 (3.80) | -3.241.1 (3.50) | -3.4113 (3.40) | -3.3114 (2.70) 
PyratBay | -1.9'1l (3.60) | -1.915 (3.30) | -1.9*1? (3.60) | -2.2453 (3.60) 
NEMESIS | -3.6'59 (3.30) | -2.1'97 (3.00) | -1.8'93 (4.10) | -4.4433 (1.80) 
SCARLET* Suy -2.5411 -3.0155 ts 
PLATON* -3.9785 -3.141 -3.5*13 -3.4*17 

log(Xcu,) (95% upper limit) 
Phase: 0.0 0.25 0.5 0.75 
HyDRA <-5.8 | <-5.3 | < —4.6 | < —56 
PyratBay | < —5.2 | < —4.4 | < —3.5 | < —4.1 
NEMESIS < —5.9 | < —4.9 | <-49 | < -1.9 
SCARLET* | < —1.7 | < —1.5 | < —1.5 | «-12 
PLATON* < —2.5 | < —3.1 | < —4.1 | < —2.6 


The H,O volume mixing ratios and CH, upper limits inferred by each retrieval framework at each phase, in log values. H,O detection significances are shown in brackets where applicable. 
Codes marked with * assume equilibrium chemistry and cannot constrain the abundances of molecules individually. For these codes, we report the abundance of water at 100 mbar. 
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Extended Data Table 3 | Chi-squared values of equilibrium chemistry retrievals with and without opacity from water 
and methane 


X? of retrieval: all molecules / no water / no methane 
Phase: 0.0 0.25 0.5 0.75 
SCARLET | 8.1/21.0/5.4 | 20.7/44.5/12.4 | 52.0/88.4/43.1 | 17.9/53.1/5.0 
PLATON | 24.8/54.6/6.1 | 43.3/84.4/46.3 | 71.1/85.7/73.2 | 43.8/73.5/47.2 


Theretrieval setup is identical aside from the inclusion or exclusion of opacity from water or methane. The results presented here do not consider the error-inflation parameter (the retrievals 
would be free to inflate errors until a x’/N,,ints value of 1 is reached). The models were fit to the phase-resolved spectra which each had 11 data points. The PLATON retrieval included 9 
parameters. The SCARLET retrieval included 7 parameters for phases 0.25, 0.5, and 0.75, and 6 parameters for phase 0.0. 
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